/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
#if !defined(_Cauchy_Heat_Equation_FD_)
#define _Cauchy_Heat_Equation_FD_

#include <set>

#include <src/math/Operators/all.h>
#include <src/algorithm/Evolution_Algorithm/Crank_Nicolson.h>
#include <src/math/CauchyProblems/CauchyProblem_traits.h>
#include <src/math/functor/GaussianFunctor.h>
#define NUMBER_OF_PARTICLES 10

/** @ingroup	heat_equations
*   @brief Definition of the Crank Nicolson Algorithm to solve the heat equation @link heat_equations@endlink in the unit cube.
*/
template <class data, std::size_t Dim>
struct CrankNicolson_heat : public Crank_Nicolson<Delta_t<data, Dim>>
{
	virtual OS_double CFL() {
		return 1./(get_state_variable()->size()*get_state_variable()->size());
	};
	virtual bool criteria() {
		OS_double breakpoint1 = jmmath::norm22(Op_iterate - Op_S_k1);
		OS_double breakpoint2 = jmmath::norm22(Op_iterate);
		return breakpoint1/(1.+breakpoint2) < 1e-4;
	};
};


/** @ingroup	heat_equations
*   @brief a finite difference scheme to solve the heat equation @link heat_equations@endlink in the unit cube.
*
*	This class initialize the Cauchy problem for the multi-dimensional heat equation
*	\f[
*		\partial_t u(t,x) = \Delta u(t,x)
*	\f];
*	with \f$x \in [0,1]^D\f$, \f$t\ge 0\f$.
*
*	It illustrate the behavior of the finite difference scheme
*	\f[
*		\delta_t u^n = \bar{\nabla} \cdot \nabla u^{n+1/2} = 0
*	\f];
*	The algorithm used is Crank Nicolson.
*	The boundary conditions are given by <data> : either Dirichlet or Periodic ones.
*	The initial data corresponds to a normal distribution centered in 0.5.
*/
template <class data, std::size_t Dim, class Interface = Graphical_Interface>
		class Cauchy_Delta : public CauchyProblem_traits< CrankNicolson_heat < data,Dim >, Interface >
		{
		public :
			virtual ~Cauchy_Delta(){}
			Cauchy_Delta(){
				initial_conditions(NUMBER_OF_PARTICLES);
			};
			virtual smart_ptr<interface_type> Create_Instance() {
				return smart_ptr<interface_type>(new Cauchy_Delta);
			};

			virtual void set_dim(OS_size dim) {
				OS_DYNAMIC_CHECK(dim == Dim,"Cauchy_Delta: only dimension Dim");
			};
			OS_size get_dim() {return Dim;};
			virtual void initial_conditions(OS_size nb_part)
			{
				std::vector<std::size_t> size_(Dim);
				std::size_t inf = std::pow(nb_part,1./Dim);
				std::fill(size_.begin(),size_.end(),inf+1);
				std::vector<double> average(Dim);
				std::fill(average.begin(),average.end(),.5);
				set_initial_conditions(Gaussian_Functor<result_type>(average)(size_));
			};
			virtual typename interface_type::p_array_2 get_representation() {
				return representation_3<data,2,Finite_Diff_Operator>()(get_solution());
			};
		};
/** @ingroup	heat_equations
*   @brief a finite difference scheme specializing equation @link heat_equations@endlink in the unit cube.
*
*	This class initialize the Cauchy problem for the one dimensional heat equation
*	\f[
*		\partial_t u(t,x) = \partial_{xx} u(t,x)
*	\f];
*	with \f$x \in [0,1]\f$, \f$t\ge 0\f$.
*
*	It illustrate the behavior of the finite difference scheme
*	\f[
*		\delta_t u^n = \delta_{x\bar{x}} u^{n+1/2} = 0
*	\f];
*	The algorithm used is Crank Nicolson.
*	The boundary conditions are given by <data> : either Dirichlet or Periodic ones.
*	The initial data corresponds to a normal distribution centered in 0.5.
*/
template <class data>
		class Cauchy_Delta<data,1,Graphical_Interface> : public CauchyProblem_traits< CrankNicolson_heat < data,1 >, Graphical_Interface >
		{
		public :
			virtual ~Cauchy_Delta(){}
			Cauchy_Delta(){
				initial_conditions(NUMBER_OF_PARTICLES);
			};
			virtual smart_ptr<interface_type> Create_Instance() {
				return smart_ptr<interface_type>(new Cauchy_Delta);
			};

			virtual void set_dim(OS_size dim) {
				OS_DYNAMIC_CHECK(dim == 1,"Cauchy_Delta: only dimension 1");
			};
			OS_size get_dim() {return 1;};
			virtual void initial_conditions(OS_size nb_part)
			{
				set_initial_conditions(Gaussian_Functor<result_type>()(nb_part));
			};
		};


#undef NUMBER_OF_PARTICLES

#endif

